Setup

library(SeqArray)
library(SNPRelate)
library(pander)
library(scales)
library(magrittr)
library(tidyverse)
library(readxl)
library(sp)
library(ggmap)
library(rgdal)
library(ggsn)
library(parallel)
library(qqman)
library(ggrepel)
library(plyranges)
library(rtracklayer)
library(AnnotationDbi)
library(GO.db)
theme_set(theme_bw())
panderOptions("missing", "")
mc <- min(12, detectCores() - 1)
alpha <- 0.05
maxKb <- 40
w <- 169/25.4
h <- 169/25.4

This workflow immediately follows 03_SNPFiltering and performs an analysis on the SNPs retained after filtering. Analysis performed is:

gdsPath <- file.path("..", "5_stacks", "gds", "populations.snps.gds")
gdsFile <- seqOpen(gdsPath, readonly = TRUE)
keepSNPs <- readRDS("keepSNPsAfterLDPruning.RDS")
sampleID <- tibble(
    Sample = seqGetData(gdsFile, "sample.annotation/Sample"),
    Population = seqGetData(gdsFile, "sample.annotation/Population"),
    Location = seqGetData(gdsFile, "sample.annotation/Location")
)
popSizes <- sampleID %>% 
    group_by(Population) %>%
    summarise(n = dplyr::n())
snps <- tibble(
    variant.id = seqGetData(gdsFile, "variant.id") ,
    chromosome = seqGetData(gdsFile, "chromosome"),
    position = seqGetData(gdsFile, "position"),
    snpID = seqGetData(gdsFile, "annotation/id")
) %>%
    mutate(
      snpID = str_remove(snpID, ":.$"),
      snpID = str_replace(snpID, ":", "_")
    ) %>%
    separate(snpID, into = c("Locus ID", "Col"), remove = FALSE) %>%
    mutate_at(vars(`Locus ID`, Col), as.integer) %>%
    mutate(
      Col = Col - 1,
      snpID = paste(`Locus ID`, Col, sep = "_")
    ) %>%
  right_join(keepSNPs)
seqSetFilter(gdsFile, variant.id = snps$variant.id)
## # of selected variants: 18,886
genotypes <- snps %>%
    cbind(
      seqGetData(gdsFile, "genotype") %>% 
        apply(MARGIN = 3, colSums) %>%
        t %>%
        set_colnames(sampleID$Sample)) %>%
  as_tibble() %>%
  gather(Sample, Genotype, one_of(sampleID$Sample)) %>%
  dplyr::filter(!is.na(Genotype)) %>%
  arrange(variant.id, Sample) %>%
  left_join(sampleID)
seqResetFilter(gdsFile)
## # of selected samples: 146
## # of selected variants: 146,763
snpIn1996 <- genotypes %>% 
    filter(Population == 1996) %>% 
    group_by(variant.id) %>% 
    summarise(maf = mean(Genotype) / 2) %>%
    filter(maf > 0)
genotypes %<>% 
    right_join(snpIn1996)
seqClose(gdsFile)

In addition to the above, the set of genes corresponding the Ensembl Release 96 were loaded.

ensGenes <- file.path("..", "external", "Oryctolagus_cuniculus.OryCun2.0.96.gff3.gz") %>%
    import.gff3(feature.type = "gene", sequenceRegionsAsSeqinfo = TRUE) %>%
    .[,c("gene_id", "Name", "description")]

The set of SNPs under investigation was also defined as a GRanges object.

snpsGR <- makeGRangesFromDataFrame(
    df = snps, 
    keep.extra.columns = TRUE, 
    ignore.strand = TRUE, 
    seqinfo = seqinfo(ensGenes), 
    seqnames.field = "chromosome", 
    start.field = "position", 
    end.field = "position"
)

This was intersected with the set of genes to connect SNPs to genes within 40 kb.

snp2Gene <- snpsGR %>% 
    resize(width = 1000*maxKb/2 - 1) %>%
    trim() %>%
    join_overlap_inner(ensGenes)

Of the initial 18,886 SNPs, 8,690 were found to be within 40KB of a gene.

A set of SNPs within exonic regions was also defined.

ensExons <- file.path("..", "external", "Oryctolagus_cuniculus.OryCun2.0.96.gff3.gz") %>%
    import.gff3(feature.type = "exon", sequenceRegionsAsSeqinfo = TRUE) %>%
  join_overlap_inner_within_directed(
    ensGenes, maxgap = 0, minoverlap = 0, suffix = c(".exon", "")
  ) %>%
  join_overlap_inner(snpsGR)
ensExons <- ensExons[,setdiff(colnames(mcols(ensExons)), "Name.exon")]

PCA

lowCall <- c("gc2901", "gc2776", "gc2731", "gc2727", "gc2686")
snp4PCA <- genotypes %>% 
    filter(!Sample %in% lowCall) %>%
    group_by(variant.id, Population) %>% 
    summarise(n = dplyr::n()) %>% 
    spread(Population, n) %>%
    mutate(N = `1996` + `2010` + `2012`) %>% 
    ungroup() %>%
    filter(N > 0.95*(sum(popSizes$n) - length(lowCall)))
pca <- genotypes %>%
    filter(variant.id %in% snp4PCA$variant.id) %>%
    dplyr::select(variant.id, Sample, Genotype) %>%
    spread(Sample, Genotype) %>%
    as.data.frame() %>%
    column_to_rownames("variant.id") %>%
    as.matrix() %>%
    apply(2, function(x){
        x[is.na(x)] <- mean(x, na.rm = TRUE)
        x
    }) %>%
    t() %>%
    .[, apply(., 2, function(x){length(unique(x)) > 1})] %>%
    prcomp( center = TRUE)

As noted in the previous section sample samples gc2901, gc2776, gc2731, gc2727 and gc2686 had a SNP identification rate \(< 50\)% and as such these were marked as potential outliers. Ignoring these samples, and restricting data to SNPs identified in \(>95\)% of all samples, a preliminary PCA was performed This amounted to 7,767 of the possible 18,878 SNPs for analysis using PCA. Missing values were specified as the mean MAF over all populations combined.

Given the initially observed structure, in which samples from the 2012 population are separating from the other samples which group with the 1996 population, the collection points for the 2012 samples as investigated.

sampleID %<>%
    left_join(
        file.path("..", "external", "GPS_Locations.xlsx") %>%
            read_excel() %>%
            dplyr::select(Sample, ends_with("tude")) %>%
            mutate(Sample = gsub("[Oo][Rr][Aa] ([0-9ABC]*)", "ora\\1", Sample))
    )
*PCA showing population structure. Point size reflects the proportion of SNPs for which imputation was required, and the observed structure appeared independent of this.*

PCA showing population structure. Point size reflects the proportion of SNPs for which imputation was required, and the observed structure appeared independent of this. 

loc <- c(
  range(sampleID$Longitude, na.rm = TRUE) %>% mean,
  range(sampleID$Latitude, na.rm = TRUE) %>% mean
)
saPoly <- readRDS(file.path("..", "external", "saPoly.RDS"))
roads <- readRDS(file.path("..", "external", "roads.RDS"))
gc <- SpatialPoints(cbind(x = 138.655972, y = -31.200305))
proj4string(gc) <- "+proj=longlat +ellps=GRS80 +no_defs"
xBreaks <- seq(138.65, 138.8, by = 0.05)
xLabs <- parse(text = paste(xBreaks, "*degree ~ E", sep = ""))
yBreaks <- seq(-31.2, -31.32, by = -0.04)
yLabs <- parse(text = paste(-yBreaks, "*degree ~ S", sep = ""))
leftN <- tibble(
  x = c(138.7965, 138.8, 138.8) - 0.01,
  y = c(-31.2, -31.198, -31.193)
)
rightN <- tibble(
  x = c(138.8, 138.8, 138.8035) - 0.01,
  y = c(-31.193, -31.198, -31.2)
)
ggMap <- get_map(loc, zoom = 12, maptype = "terrain", color = "bw")
zoomPlot <- ggmap(ggMap, extent = "normal", maprange = FALSE) +
    geom_path(
        aes(long, lat, group = group), 
        data = subset(roads, SURFACE == "UNSE"), 
        linetype = 2, size = 0.3) + 
    geom_path(
        aes(long, lat, group = group), 
        data = subset(roads, SURFACE != "UNSE"), 
        linetype = 1, size = 0.4) +
    geom_label(x = 138.74, y = -31.29, label = "Flinders Ranges NP", alpha = 0.4) +
    geom_label(x = 138.72, y = -31.22, label = "Gum Creek", alpha = 0.4) +
    geom_point(aes(x, y), data = as.data.frame(gc), shape = 3, size = 3) +
    geom_text(aes(x, y), label = "HS", data = as.data.frame(gc), nudge_y = 0.005) +
    geom_polygon(
        data = subset(saPoly, F_CODE == "HD"),
        aes(long, lat, group = group),
        fill = rgb(1, 1, 1, 0), colour = "grey10", size = 0.3) +
    geom_polygon(
        data = subset(saPoly, F_CODE == "PARK"),
        aes(long, lat, group = group),
        fill = rgb(1, 1, 1, 0), colour = "grey10", size = 0.2) +
  geom_point(
    data= filter(pca4Plot, grepl("2012", Population)),
    aes(Longitude, Latitude, colour = Population),
    size = 0.9*ps) +
  geom_polygon(
    aes(x, y), data = leftN, fill = "white", colour = "grey10", size = 0.4
  ) +
  geom_polygon(
    aes(x, y), data = rightN, fill = "grey10", colour = "grey10", size = 0.4
  ) +
  geom_text(
    x = 138.79, y = -31.19, label = "N", colour = "grey10", size = 4
  ) +
  scale_colour_manual(values = popCols[2:3]) +
  coord_cartesian(
    xlim = c(138.618, 138.81),
    ylim = c(-31.335, -31.18),
    expand = 0
  ) +
  scale_x_continuous(breaks = xBreaks, labels = xLabs) +
  scale_y_continuous(breaks = yBreaks, labels = yLabs) +
  guides(colour = FALSE) +
  ggsn::scalebar(
    x.min = 138.618,
    x.max = 138.81,
    y.min = -31.335, 
    y.max = -31.18,
    transform = TRUE,
    dist = 2, 
    dist_unit = "km",
    model = 'GRS80',
    height = 0.012, 
    st.size = 4,
    location = 'bottomright',
    anchor = c(x = 138.8, y = -31.328)) +
  labs(x = "Longitude", y = "Latitude") +
  theme(
    text = element_text(size = fs),
    plot.margin = unit(c(1, 1, 1, 1), "mm")
  )
ausPolygon <- readRDS(file.path("..", "external", "ausPolygon.RDS"))
ausPts <- SpatialPoints(cbind(x = loc[1], y = loc[2]))
proj4string(ausPts) <- proj4string(ausPolygon)
ausPlot <- ggplot() +
  geom_polygon(
    data = ausPolygon, 
    aes(long, lat, group = group),
    fill = "white", colour = "black"
  ) + 
  geom_point(data = as.data.frame(ausPts), aes(x, y), size = 1.5) +
  theme_void() +
  theme(plot.background = element_rect(fill = "white", colour = "black"))
## png 
##   2
*Figure 1: Collection points for all 2012 samples with colours showing sub-populations initially defined by PCA analysis and *k*-means clustering.*

Figure 1: Collection points for all 2012 samples with colours showing sub-populations initially defined by PCA analysis and k-means clustering. 

zoomLoc <- c(138.753, -31.242)
xBreaks <- seq(138.74, 138.76, by = 0.01)
xLabs <- parse(text = paste(xBreaks, "*degree ~ W", sep = ""))
yBreaks <- seq(-31.235, -31.25, by = -0.005)
yLabs <- parse(text = paste(-yBreaks, "*degree ~ S", sep = ""))
central <- rbind(
  x = c(138.749, 138.755),
  y = c(-31.2365, -31.2495)
) %>%
  set_colnames(c("min", "max"))
leftN <- tibble(
  x = c(138.7595, 138.76, 138.76),
  y = c(-31.234, -31.2335, -31.2325)
)
rightN <- tibble(
  x = c(138.76, 138.76, 138.7605),
  y = c(-31.2325, -31.2335, -31.234)
)
get_map(zoomLoc, zoom = 15, maptype = "terrain", source = "google", color = "bw") %>%
  ggmap() +
  geom_jitter(
    data = filter(pca4Plot, grepl("2012", Population)), 
    aes(Longitude, Latitude, colour = Population), 
    size = 3, width = 0.0005, height = 0) +
  geom_rect(
    xmin = central["x", "min"],
    xmax = central["x", "max"],
    ymin = central["y", "min"],
    ymax = central["y", "max"],
    fill = "red", alpha = 0.01, colour = "black") +
  geom_polygon(
    aes(x, y), data = leftN, fill = "white", colour = "grey10", size = 0.4
  ) +
  geom_polygon(
    aes(x, y), data = rightN, fill = "grey10", colour = "grey10", size = 0.4
  ) +
  geom_text(
    x = 138.76, y = -31.232, label = "N", colour = "grey10", size = 5
  ) +
  scale_colour_manual(values = popCols[2:3]) +
  scale_x_continuous(breaks = xBreaks, labels = xLabs) +
  scale_y_continuous(breaks = yBreaks, labels = yLabs) +
  theme_bw() +
  guides(colour = FALSE) +
  labs(x = "Longitude", y = "Latitude") +
  coord_cartesian(
    xlim = c(138.74, 138.762),
    ylim = c(-31.253, -31.231),
    expand = 0
  ) +
  ggsn::scalebar(
    x.min = 138.74, x.max = 138.762, 
    y.min = -31.253, y.max = -31.231,
    transform = TRUE,
    dist = 0.25, dist_unit = "km",, model = 'GRS80', 
    height = 0.012, st.size = 4,
    location = 'bottomright',
    anchor = c(x = 138.761, y = -31.252)
    )
*Zoomed-in view of the central region for 2012 samples with colours showing sub-populations defined by PCA analysis. The region considered to be the Central Region is shaded in red. Due to overlapping GPS points a small amount of jitter has been added to the x-axis.*

Zoomed-in view of the central region for 2012 samples with colours showing sub-populations defined by PCA analysis. The region considered to be the Central Region is shaded in red. Due to overlapping GPS points a small amount of jitter has been added to the x-axis.

Region Analysis

Removal of SNPs Associated with Collection Region

The structure observed within the 2012 population in the PCA could possibly be explained by recent migration into this region. As the samples collected in the outer regions appeared very similar to the 1996 population in the above plots, this would possibly indicate migration a very recent event as the genetic influence of this has not spread through the wider area. Although this may be due to other factors such as sampling bias, this structure was addressed by identifying SNPs which showed an association with the sub-populations identified by PCA analysis. In this way, any candidate SNPs obtained below will be less impacted by this structure, and will be more reflective of the intended variable under study, i.e. selection over time, as opposed to any internal structure of the 2012 population.

oraRegions <- pca4Plot %>% 
  filter(grepl("2012", Population)) %>% 
  rowwise() %>%
  mutate(
    yCentral = cut(Latitude, breaks = central["y",], include.lowest = TRUE), 
    xCentral = cut(Longitude, breaks = central["x",], include.lowest = TRUE),
    Central = (is.na(yCentral) + is.na(xCentral)) == 0
  ) %>%
  dplyr::select(Sample, Central)

Testing for Structure in 2012

This model tests:
H0: No association between genotypes and collection region
HA: An association exists between genotypes and collection region

regionResults <- genotypes %>%
    filter(Population == 2012) %>%
    split(f = .$variant.id) %>%
    mclapply(
      function(x){
        ft <- list(p.value = NA)
        if (length(unique(x$Genotype)) > 1) {
          ft <- x %>%
            left_join(oraRegions) %>%
            group_by(Genotype, Central) %>%
            tally() %>%
            spread(Genotype, n, fill = 0) %>%
            column_to_rownames("Central") %>%
            fisher.test()
        }
        x %>%
          distinct(variant.id, snpID, chromosome, position) %>%
          mutate(p = ft$p.value)
      }, 
      mc.cores = mc
    ) %>%
  bind_rows() %>%
  filter(!is.na(p)) %>%
  arrange(p)

A total of 1682 SNPs were detected as showing a significant association between genotype and the collection region. Under H0, the number expected using α = 0.05 would be 943, and as this number was approximately double that expected, this was taken as evidence of this being a genuine point of concern for this dataset.

Notably, Type II errors were of principle concern in this instance, and as such every SNP with p < 0.05 in the above test was excluded from downstream analysis.

regionSNPs <- filter(regionResults, p < 0.05)
saveRDS(regionSNPs, "regionSNPs.RDS")

Under this additional filtering step, the original set of 18784 SNPs will be reduced to 17,102 for testing by genotype and allele frequency.

Verification Of Removal

In order to verify that the removal of the above SNPs removed the undesired population structure from the 2012 population, the above PCA was repeated, excluding the SNPs marked for removal. The previous structure noted in the data was no longer evident, and as such, these SNPs were marked for removal during analysis by genotype and allele frequency.

pcaPost <- genotypes %>%
  filter(
    variant.id %in% snp4PCA$variant.id,
    !variant.id %in% regionSNPs$variant.id
  ) %>%
  dplyr::select(variant.id, Sample, Genotype) %>%
  spread(Sample, Genotype) %>%
  as.data.frame() %>%
  column_to_rownames("variant.id") %>%
  as.matrix() %>%
  apply(2, function(x){
    x[is.na(x)] <- mean(x, na.rm = TRUE)
    x
  }) %>%
  t() %>%
  .[, apply(., 2, function(x){length(unique(x)) > 1})] %>%
  prcomp( center = TRUE) 
pcaPost4Plot <- pcaPost$x %>%
  as.data.frame() %>%
  rownames_to_column("Sample") %>%
  as_tibble() %>%
  dplyr::select(Sample, PC1, PC2, PC3) %>%
  left_join(sampleID) %>%
  mutate(Cluster = kmeans(cbind(PC1, PC2, PC3), 3)$cluster) %>%
  group_by(Cluster) %>%
  mutate(maxY = max(Latitude, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(
    Population = case_when(
      Population == 1996 ~ "1996",
      Population == 2010 ~ "Outgroup (Turretfield)",
      maxY == max(maxY) ~ "2012 (Outer)",
      maxY != max(maxY) ~ "2012 (Central)"
    )
  ) %>%
  left_join(
    genotypes %>% 
      filter(variant.id %in% snp4PCA$variant.id, !Sample %in% lowCall) %>% 
      group_by(Sample) %>% 
      tally() %>% 
      mutate(imputationRate = 1 - n / nrow(snp4PCA))
  )
*Figure 2: Principal Components Analysis showing structures before removal of SNPs denoting collection region in the 2012 population, and after removal of these SNPs*

Figure 2: Principal Components Analysis showing structures before removal of SNPs denoting collection region in the 2012 population, and after removal of these SNPs

SNP Analysis

Genotype Frequency Model

This model tests:
H0: No association between genotypes and populations
HA: An association exists between genotypes and populations

genotypeResults <- genotypes %>%
  filter(
    Population != 2010,
    !variant.id %in% regionSNPs$variant.id,
    variant.id %in% snp2Gene$variant.id
  ) %>%
  group_by(variant.id, snpID, Population, Genotype) %>%
  tally() %>%
  ungroup() %>%
  split(f = .$variant.id) %>%
  mclapply(function(x){
    ft <- list(p.value = NA)
    if (length(unique(x$Genotype)) > 1) {
      ft <- x %>%
        spread(Genotype, n, fill = 0) %>%
        column_to_rownames("Population") %>%
        dplyr::select(-variant.id, -snpID) %>%
        fisher.test()
    }
    x %>% 
      distinct(variant.id) %>%
      mutate(p = ft$p.value)
  },
  mc.cores = mc
  ) %>%
  bind_rows() %>%
  filter(!is.na(p)) %>%
  mutate(FDR = p.adjust(p, "fdr"), adjP = p.adjust(p, "bonferroni")) %>%
  arrange(p) %>%
  left_join(
    genotypes %>%
      distinct(variant.id, snpID, chromosome, position)
  ) %>%
  dplyr::select(variant.id, snpID, chromosome, position, everything())
genoSNPs <- genotypeResults %>%
  dplyr::filter(FDR < alpha) %>%
  .[["snpID"]]
Candidate SNPs to an FDR of 5% when analysing by genotype. The nearest gene to each SNP (within 40kb ) is indicated. Bonferroni-adjusted p-values are also given in the final column.
snpID chromosome position gene_id Name dist2Gene p FDR adjP
686773_29 3 90,851,512 ENSOCUG00000009616 CRISPLD1 CDS 1.053e-07 0.0006831 0.0008341
916731_91 4 89,602,973 ENSOCUG00000006126 RIC8B Intronic 1.865e-07 0.0006831 0.001478
836950_151 4 35,111,855 ENSOCUG00000012144 TMPRSS12 5,745 2.587e-07 0.0006831 0.002049
2227006_109 13 129,650,867 ENSOCUG00000022823 ID3 16,426 4.419e-07 0.0007665 0.0035
3040789_114 19 30,689,426 ENSOCUG00000010440 MSI2 Intronic 4.838e-07 0.0007665 0.003832
3695415_108 GL018751 737,318 ENSOCUG00000003672 TRAF3 Intronic 1.837e-06 0.002369 0.01455
1902686_98 12 60,113,885 ENSOCUG00000001871 12,820 2.094e-06 0.002369 0.01658
644062_1 3 56,309,981 ENSOCUG00000004827 SFXN1 Intronic 3.509e-06 0.003475 0.0278
4149102_107 GL019077 34,265 ENSOCUG00000024268 IFT140 Intronic 4.816e-06 0.003913 0.03814
87917_119 1 64,635,286 ENSOCUG00000007758 CEP78 Intronic 5.487e-06 0.003913 0.04346
2689075_49 16 56,527,308 ENSOCUG00000000984 ESRRG Intronic 5.622e-06 0.003913 0.04454
4010189_98 GL018907 283,766 ENSOCUG00000008870 BRWD1 Intronic 5.928e-06 0.003913 0.04696
4044777_119 GL018933 204,058 ENSOCUG00000005859 FNBP1 Intronic 7.754e-06 0.004725 0.06142
4146664_90 GL019084 74,038 ENSOCUG00000010071 GNB1 Intronic 1.842e-05 0.01042 0.1459
3016258_99 19 19,178,148 ENSOCUG00000015254 SEZ6 Intronic 2.168e-05 0.0108 0.1717
4176850_23 GL019154 883 ENSOCUG00000002484 MTHFSD 828 2.181e-05 0.0108 0.1728
4098854_101 GL018985 129,495 ENSOCUG00000006396 Intronic 3.086e-05 0.01438 0.2444
3829050_109 GL018791 957,960 ENSOCUG00000007849 FHAD1 Intronic 4.006e-05 0.01704 0.3173
482617_100 2 139,723,776 ENSOCUG00000001127 RHOQ 8,346 4.335e-05 0.01704 0.3434
3378006_1 GL018704 2,271,135 ENSOCUG00000006255 AGO3 9,610 4.449e-05 0.01704 0.3524
2580149_42 15 87,473,844 ENSOCUG00000001910 ADGRL3 Intronic 4.518e-05 0.01704 0.3579
4136275_162 GL019049 178,272 ENSOCUG00000026644 Intronic 9.113e-05 0.03281 0.7219
2162055_71 13 91,091,107 ENSOCUG00000003017 ERICH3 Intronic 0.0001249 0.04301 0.9892
1287228_44 8 6,997,901 ENSOCUG00000017747 CPNE8 Intronic 0.0001431 0.04542 1
3752363_81 GL018765 531,530 ENSOCUG00000022152 RFC2 12,675 0.0001434 0.04542 1
870312_62 4 57,287,528 ENSOCUG00000004936 NAV3 Intronic 0.0001566 0.04772 1
3812992_139 GL018789 372,825 ENSOCUG00000008689 Intronic 0.0001653 0.04849 1

Under the full genotype model:

  • 12 genotypes were detected as being significantly associated with the two populations when controlling the FWER at α = 0.05
  • 27 genotypes were detected as being significantly associated with the two populations when controlling the FDR at α = 0.05
  • For the most highly ranked SNP (686773_29), the minor allele has been completely lost in the 2012 population
Figure 3a: Manhattan plot showing results for all SNPs on chromosomes 1:21 for analysis by genotype. The horizontal line indicates the cutoff for an FDR of 5%, with SNPs with an adjusted p-value < 0,05 shown in green

Figure 3a: Manhattan plot showing results for all SNPs on chromosomes 1:21 for analysis by genotype. The horizontal line indicates the cutoff for an FDR of 5%, with SNPs with an adjusted p-value < 0,05 shown in green

*Variants detected as significant to an FDR of 5%. Those denoted with an asterisk received a Bonferroni-adjusted p-value < 0.05 and these typically involved reduction of the minor allele frequency and a shift towards homozygous reference. The remaining sites showed a combination of the same and increasing minor allele abundance, with some variants becoming exclusively heterozygous in the 2012 population.*

Variants detected as significant to an FDR of 5%. Those denoted with an asterisk received a Bonferroni-adjusted p-value < 0.05 and these typically involved reduction of the minor allele frequency and a shift towards homozygous reference. The remaining sites showed a combination of the same and increasing minor allele abundance, with some variants becoming exclusively heterozygous in the 2012 population. 

write_tsv(genotypeResults, "genotypeResults.tsv")

Allele Frequency Model

This model tests:
H0: No association between allele frequencies and populations
HA: An association exists between allele frequencies and populations

alleleResults <- genotypes %>%
  filter(
    Population != 2010,
    !variant.id %in% regionSNPs$variant.id,
    variant.id %in% snp2Gene$variant.id
  ) %>%
  group_by(snpID, Population) %>%
  summarise(
    P = sum(2 - Genotype),
    Q = sum(Genotype)
  ) %>%
  ungroup() %>%
  split(f = .$snpID) %>%
  mclapply(function(x){
    m <- as.matrix(x[c("P", "Q")])
    ft <- list(p.value = NA) 
    if (length(m) == 4) {
      ft <- fisher.test(m)
    }
    x %>%
      mutate(MAF = Q / (P + Q)) %>%
      dplyr::select(snpID, Population, MAF) %>%
      spread(Population, MAF) %>%
      mutate(p = ft$p.value)
  },mc.cores = mc) %>%
  bind_rows() %>%
  filter(!is.na(p)) %>%
  dplyr::rename(
    MAF_1996 = `1996`,
    MAF_2012 = `2012`
  ) %>%
  mutate(
    FDR = p.adjust(p, "fdr"),
    adjP = p.adjust(p, "bonferroni"),
    # Make the MAF is relative to the 1996 population not the reference
    MAF_2012 = case_when(
      MAF_1996 > 0.5 ~ 1 - MAF_2012,
      MAF_1996 <= 0.5 ~ MAF_2012
    ),
    MAF_1996 = case_when(
      MAF_1996 > 0.5 ~ 1 - MAF_1996,
      MAF_1996 <= 0.5 ~ MAF_1996
    )
  ) %>%
  arrange(p) %>%
  left_join(
    genotypes %>% distinct(snpID, chromosome, position)
  ) %>%
  dplyr::select(snpID, chromosome, position, everything())
alleleSnps <- alleleResults %>%
  dplyr::filter(FDR < alpha) %>%
  .[["snpID"]]
SNPs considered as significant when analysing by genotype using an FDR cutoff of 5% Minor allele frequencies are reported based on frequencies in the 1996 population. Any genes within 40kb are also shown. Bonferroni-adjusted p-values are given in the final column.
snpID chromosome position MAF_1996 MAF_2012 gene_id Name dist2Gene p FDR adjP
916731_91 4 89,602,973 0.1961 0 ENSOCUG00000006126 RIC8B Intronic 6.875e-07 0.005452 0.005452
3040789_114 19 30,689,426 0.1863 0 ENSOCUG00000010440 MSI2 Intronic 1.531e-06 0.006068 0.01214
836950_151 4 35,111,855 0.2907 0.04082 ENSOCUG00000012144 TMPRSS12 5,745 3.086e-06 0.007532 0.02447
2227006_109 13 129,650,867 0.2708 0.03333 ENSOCUG00000022823 ID3 16,426 4.172e-06 0.007532 0.03308
1902686_98 12 60,113,885 0.1633 0 ENSOCUG00000001871 12,820 4.749e-06 0.007532 0.03766
4044777_119 GL018933 204,058 0.01754 0.2019 ENSOCUG00000005859 FNBP1 Intronic 7.749e-06 0.01024 0.06145
3695415_108 GL018751 737,318 0.2609 0.03261 ENSOCUG00000003672 TRAF3 Intronic 1.275e-05 0.0131 0.1011
87917_119 1 64,635,286 0.186 0 ENSOCUG00000007758 CEP78 Intronic 1.43e-05 0.0131 0.1134
4010189_98 GL018907 283,766 0.1932 0.01 ENSOCUG00000008870 BRWD1 Intronic 1.627e-05 0.0131 0.129
2689075_49 16 56,527,308 0.05814 0.3043 ENSOCUG00000000984 ESRRG Intronic 1.652e-05 0.0131 0.131
4149102_107 GL019077 34,265 0.2232 0.02941 ENSOCUG00000024268 IFT140 Intronic 2.279e-05 0.01609 0.1807
1997286_93 12 141,844,290 0.08889 0.3556 ENSOCUG00000004829 ESR1 Intronic 2.435e-05 0.01609 0.1931
2405359_66 14 97,813,263 0.02083 0.2021 ENSOCUG00000005835 STXBP5L Intronic 4.48e-05 0.02733 0.3553
4146664_90 GL019084 74,038 0.1961 0.02041 ENSOCUG00000010071 GNB1 Intronic 5.49e-05 0.0311 0.4354
4098854_101 GL018985 129,495 0.1702 0.01042 ENSOCUG00000006396 Intronic 6.686e-05 0.03214 0.5302
4098838_14 GL018985 123,328 0.3636 0.1275 ENSOCUG00000006396 Intronic 6.797e-05 0.03214 0.539
1045185_75 7 2,702,919 0.2143 0.03774 ENSOCUG00000027861 3,276 7.807e-05 0.03214 0.6191
3016258_99 19 19,178,148 0.22 0.03261 ENSOCUG00000015254 SEZ6 Intronic 7.861e-05 0.03214 0.6234
3829050_109 GL018791 957,960 0.163 0.009615 ENSOCUG00000007849 FHAD1 Intronic 8.094e-05 0.03214 0.6419
2580149_42 15 87,473,844 0.3511 0.1154 ENSOCUG00000001910 ADGRL3 Intronic 8.416e-05 0.03214 0.6674
3999128_75 GL018883 283,912 0.3654 0.125 ENSOCUG00000016318 SERPINA6 Intronic 8.683e-05 0.03214 0.6885
4176850_23 GL019154 883 0.234 0.04082 ENSOCUG00000002484 MTHFSD 828 8.916e-05 0.03214 0.707
686773_29 3 90,851,512 0.4375 0.1667 ENSOCUG00000009616 CRISPLD1 CDS 9.731e-05 0.03355 0.7716
1997906_1 12 142,355,344 0.1556 0.413 ENSOCUG00000023878 SYNE1 Intronic 0.0001422 0.04699 1
3455893_29 GL018713 365,938 0.2264 0.4889 ENSOCUG00000009435 XKR5 Intronic 0.0001552 0.04922 1

Under this model:

  • 5 SNP alleles were detected as being significantly associated with the two populations when controlling the FWER at α = 0.05.
  • 25 SNP alleles were detected as being significantly associated with the two populations when controlling the FDR at α = 0.05
## png 
##   2
Manhattan plot showing results for all SNPs on chromosomes 1:21 when analysing by allele frequencies. The horizontal line indicates the cutoff for an FDR of 5%, with SNPs considered significant under the Bonferroni adjustment shown in green.

Manhattan plot showing results for all SNPs on chromosomes 1:21 when analysing by allele frequencies. The horizontal line indicates the cutoff for an FDR of 5%, with SNPs considered significant under the Bonferroni adjustment shown in green.

*Comparison of minor allele frequencies in both populations. SNPs with an FDR-adjusted p-value < 0.05 are coloured red, whilst those with a Bonferroni-adjusted p-value are additionally labelled.*

Comparison of minor allele frequencies in both populations. SNPs with an FDR-adjusted p-value < 0.05 are coloured red, whilst those with a Bonferroni-adjusted p-value are additionally labelled.

*Variants detected as significant to an FDR of 5% under the Allele Frequency model. Those denoted with an asterisk received a Bonferroni-adjusted p-value < 0.05 and these typically involved reduction of the minor allele frequency and a shift towards homozygous reference. The remaining sites showed a combination of the same and increasing minor allele abundance.*

Variants detected as significant to an FDR of 5% under the Allele Frequency model. Those denoted with an asterisk received a Bonferroni-adjusted p-value < 0.05 and these typically involved reduction of the minor allele frequency and a shift towards homozygous reference. The remaining sites showed a combination of the same and increasing minor allele abundance. 

write_tsv(alleleResults, "alleleResults.tsv")

Analysis Using the FLK model

This was performed separately and no SNPs of specific interest were detected.

Export of Data for Bayescan

A VCF was required with only the 1996 and 2012 populations, and restricted to the candidate SNPs after pruning for linkage disequilibrium and detection of the allele in the 1996 population.

var2VCF <- genotypes %>% 
  distinct(variant.id, snpID) %>%
  filter(snpID %in% filter(regionResults, p > 0.05)$snpID) %>% 
  .[["variant.id"]]
gdsFile <- seqOpen(gdsPath, readonly = TRUE)
seqSetFilter(
  gdsFile, 
  variant.id = var2VCF,
  sample.id = sampleID %>% 
    filter(Population %in% c(1996, 2012)) %>% 
    .[["Sample"]]
)
seqGDS2VCF(gdsFile, "../5_stacks/vcf/filtered.vcf.gz")
seqResetFilter(gdsFile)
seqClose(gdsFile)

Inclusion of Bayescan Results

After running Bayescan, results were imported. Unfortunately, these results had IDs which were taken from the line in the above VCF and need to be converted back to IDs compatible with previous analysis.

bayRes <- file.path("..", "external", "BayescanOut.txt.gz") %>%
  gzfile() %>%
  read_delim(delim = " ", col_names = FALSE, skip = 1, col_types = "innnnn-") %>%
  set_colnames(c("ID", "prob","log10.PO.","qval","alpha","fst")) %>%
  mutate(variant.id = var2VCF[ID])
bayRes %>%
  dplyr::select(variant.id, qval, fst) %>%
  dplyr::filter(
    qval < alpha,
    variant.id %in% snp2Gene$variant.id
  ) %>%
  left_join(
    genotypes %>% 
      distinct(snpID, .keep_all = TRUE) %>% 
      dplyr::select(variant.id, chromosome, position, snpID)
  ) %>%
  makeGRangesFromDataFrame(
    keep.extra.columns = TRUE,
    ignore.strand = TRUE,
    seqinfo = seqinfo(snp2Gene), 
    seqnames.field = "chromosome",
    start.field = "position",
    end.field = "position"
  ) %>%
  join_nearest(ensGenes) %>%
  as.data.frame() %>%
  left_join(
    ensGenes %>% as.data.frame(),
    by = c("Name", "gene_id", "seqnames", "description"),
    suffix = c("_snp", "_gene")
  ) %>%
  as_tibble() %>%
  dplyr::select(seqnames, start_snp, variant.id, snpID, gene_id, Name, start_gene, end_gene, qval, fst) %>%
  mutate(dist2Gene = case_when(
    start_snp >= start_gene & start_snp <= end_gene ~ 0L,
    start_snp < start_gene ~ start_gene - start_snp,
    start_snp > end_gene ~ start_snp - end_gene
  ),
  Name = case_when(
    dist2Gene > maxKb*1000 ~ "",
    dist2Gene < maxKb*1000 ~ Name
  ),
    gene_id = case_when(
    dist2Gene > maxKb*1000 ~ "",
    dist2Gene < maxKb*1000 ~ gene_id
  ),
  inExon = snpID %in% ensExons$snpID,
  dist2Gene = case_when(
    Name == "" ~ "",
    inExon ~ "CDS",
    dist2Gene == 0 & !inExon ~ "Intronic",
    dist2Gene > 0 ~ comma(dist2Gene)
  ),
  seqnames = as.character(seqnames)
  ) %>%
  dplyr::select(
    snpID, 
    chromosome = seqnames, position = start_snp, 
    gene_id, Name, dist2Gene, fst, qval) %>%
  mutate(position = comma(position)) %>%
  arrange(qval) %>%
  pander(
    justify = "llrllrrr",
    style = "rmarkdown",
     big.mark = ",",
    split.tables = Inf,
      caption = paste(
        "Results from Bayescan with a q-value <", paste0(alpha, "."),
        "The nearest gene within", paste0(maxKb, "kb"), "is also given.",
        "One SNP detected by Bayescan was located on a scaffold with no genes (GL019119) & is not shown in this table."
      )
  )
Results from Bayescan with a q-value < 0.05. The nearest gene within 40kb is also given. One SNP detected by Bayescan was located on a scaffold with no genes (GL019119) & is not shown in this table.
snpID chromosome position gene_id Name dist2Gene fst qval
836950_151 4 35,111,855 ENSOCUG00000012144 TMPRSS12 5,745 0.03523 0.0182
2227006_109 13 129,650,867 ENSOCUG00000022823 ID3 16,426 0.03495 0.02261
4044777_119 GL018933 204,058 ENSOCUG00000005859 FNBP1 Intronic 0.03281 0.03061
3695415_108 GL018751 737,318 ENSOCUG00000003672 TRAF3 Intronic 0.03319 0.03836
4010189_98 GL018907 283,766 ENSOCUG00000008870 BRWD1 Intronic 0.03274 0.04401
*Summary of results from all three methods, using an FDR of 0.05 for each approach. All SNPs detected under Bayescan were identified using Fisher's Exact Test for both allele and genotype analysis*

Summary of results from all three methods, using an FDR of 0.05 for each approach. All SNPs detected under Bayescan were identified using Fisher’s Exact Test for both allele and genotype analysis 

sigSnps <- c(
  filter(alleleResults, FDR < alpha)$snpID,
  filter(genotypeResults, FDR < alpha)$snpID
) %>%
  unique()

Export of values for Minotaur

In order to run Minotaur a set of values was required for each SNP and a file containing the following values was prepared:

  • p-values from both the genotype & allele frequency tests above
  • the values prob, qval, alpha and fst from the Bayescan results. Only one of these should be submitted to Minotaur.
  • Overall Pi from the stacks populations output.
genotypes %>% 
  distinct(variant.id, chromosome, position, snpID) %>% 
  filter(
    !snpID %in% regionSNPs$snpID,
    variant.id %in% snp2Gene$variant.id
  ) %>%
  left_join(
    genotypeResults %>%
      dplyr::select(variant.id, chromosome, position, snpID, genotype_p = p)
  ) %>%
  left_join(
    alleleResults %>%
      dplyr::select(chromosome, position, snpID, allele_p = p)
  ) %>%
  left_join(
    bayRes %>%
      dplyr::select(
        variant.id, 
        bayescan_prob = prob, 
        bayescan_qval = qval, 
        bayescan_alpha = alpha, 
        bayescan_fst = fst)
  ) %>%
  left_join(
    file.path("..", "5_stacks", "stacks", "populations.fst_1-2.tsv.gz") %>%
  gzfile() %>%
  read_tsv(col_types = "c--ciinnnnnnnnnnn") %>%
  mutate(snpID = paste(`# Locus ID`, Column, sep = "_")) %>%
  dplyr::select(
    chromosome = Chr,
    position = BP,
    snpID,
    populations_pi = `Overall Pi`
  )
  ) %>%
  write_tsv("valuesForMinotaur.tsv")

Enrichment testing

fullPaths <- list(
  toTable(GOBPANCESTOR) %>% set_names(c("go_id", "ancestor")),
  toTable(GOCCANCESTOR) %>% set_names(c("go_id", "ancestor")),
  toTable(GOMFANCESTOR) %>% set_names(c("go_id", "ancestor"))
) %>% 
  bind_rows() %>%
  dplyr::filter(!ancestor %in% c("all", "GO:0008150", "GO:0003674", "GO:0005575")) %>%
  as_tibble() %>%
  bind_rows(
    tibble(
      go_id = unique(.$go_id),
      ancestor = unique(.$go_id)
    )
  ) %>%
  arrange(go_id, ancestor)
goSummaries <- url("https://uofabioinformaticshub.github.io/summaries2GO/data/goSummaries.RDS") %>%
  readRDS() %>%
  dplyr::rename(go_id = id)

GO Analysis By Gene

dist <- 1000*maxKb

For the initial GO analysis, any gene within 40kb of a SNP considered as significant were tested for any enrichment of biological characteristics.

allSnpsGR <- snps %>%
  dplyr::select(chromosome, position, snpID) %>%
  makeGRangesFromDataFrame(
    keep.extra.columns = TRUE,
    seqnames.field = "chromosome", 
    start.field = "position", 
    end.field = "position", 
    seqinfo = seqinfo(ensGenes)) %>%
  join_overlap_inner(
    ensGenes,
    maxgap = dist
  )
sigSnpGR <- allSnpsGR %>%
  subset(snpID %in% sigSnps)
nGenes <- allSnpsGR$gene_id %>%
  unique() %>%
  length()
nSigGenes <- sigSnpGR$gene_id %>%
  unique() %>%
  length()

10821 genes were identified as being within 40kb of a SNP. For the 34 candidate SNPs within this distance of a gene, 58 genes were identified, and these were analysed for enrichment of any GO terms.

Mappings from genes to GO terms were downloaded manually from the biomart server at www.ensembl.org.

bm <- file.path("..", "external", "gene2GO_biomart.tsv.gz") %>%
  gzfile() %>%
  read_tsv() %>%
  dplyr::filter(!is.na(`GO domain`)) %>%
  dplyr::rename(gene_id = `Gene stable ID`,
         ontology = `GO domain`,
         go_id = `GO term accession`,
         go_term = `GO term name`) %>%
  dplyr::filter(gene_id %in% allSnpsGR$gene_id) %>%
  left_join(fullPaths) %>%
  dplyr::select(gene_id, go_id = ancestor) %>%
  distinct(gene_id, go_id) %>%
  dplyr::filter(!is.na(go_id))

Complete paths of GO terms back to ontology roots were obtained for each gene. The minimum number of steps back to the ontology roots were also obtained for each GO ID.

goRes <- bm %>%
  mutate(inSig = gene_id %in% sigSnpGR$gene_id) %>%
  group_by(go_id) %>%
  summarise(N = dplyr::n(), 
            nSig = sum(inSig)) %>%
  left_join(goSummaries) %>%
  filter(nSig > 0, shortest_path > 2) %>%
  distinct(go_id, .keep_all = TRUE) %>%
  split(f = .$go_id) %>%
  mclapply(function(x){
    ft <- matrix(
      c(nSigGenes - x$nSig, x$nSig, nGenes - nSigGenes - x$N, x$N),
      nrow = 2, byrow = TRUE) %>%
      fisher.test()
    mutate(x, p = ft$p.value)
  }, mc.cores = mc) %>%
  bind_rows() %>%
  mutate(Expected = nSigGenes * N / nGenes,
         FDR = p.adjust(p, "fdr"),
         adjP = p.adjust(p, "bonferroni"),
         term = Term(go_id),
         ontology = Ontology(go_id)) %>%
  dplyr::select(go_id, term, ontology, N, nSig, Expected, p, adjP, FDR) %>%
  arrange(p)
goRes %>%
  dplyr::filter(FDR == min(FDR), nSig > 1) %>%
  mutate(Term = Term(go_id),
         Genes  = vapply(go_id, function(x){
           subset(sigSnpGR, gene_id %in% I(
             dplyr::filter(bm, go_id == x) %>%
               .[["gene_id"]]))$Name %>%
             sort() %>%
             unique() %>%
             paste(collapse = ", ")
         }, character(1))) %>%
  dplyr::select(ID = go_id, Ontology = ontology, Term, Total = N, Sig = nSig, Expected, p, FDR, Genes) %>%
  mutate(ID = gsub(":", "\\\\:", ID),
         Expected = round(Expected, 3),
         p = round(p, 4),
         FDR = round(FDR, 3)) %>%
  distinct(Ontology, Genes, .keep_all = TRUE) %>%
  pander(split.tables = Inf, style = "rmarkdown", 
         justify = "lllrrrrrl",
         caption = paste("Most highly ranked GO terms.", 
                         "Where multiple terms resulted from the identical combination of genes, only the most highly ranked terms is shown.",
                         "Overall, this corresponds to an FDR of", 
                         percent(max(.$FDR))))
Most highly ranked GO terms. Where multiple terms resulted from the identical combination of genes, only the most highly ranked terms is shown. Overall, this corresponds to an FDR of 29.5%
ID Ontology Term Total Sig Expected p FDR Genes
GO:0005496 MF steroid binding 34 3 0.182 0.001 0.295 ESR1, ESRRG, SERPINA6
GO:0002221 BP pattern recognition receptor signaling pathway 51 3 0.273 0.003 0.295 ESR1, TRAF3, UFD1
GO:1990752 CC microtubule end 17 2 0.091 0.0046 0.295 CLIP2, NAV3
GO:0003727 MF single-stranded RNA binding 19 2 0.102 0.0056 0.295 AGO3, MSI2
GO:0031050 BP dsRNA processing 19 2 0.102 0.0056 0.295 AGO3, ESR1
GO:0004879 MF nuclear receptor activity 20 2 0.107 0.0061 0.295 ESR1, ESRRG
GO:0050688 BP regulation of defense response to virus 20 2 0.107 0.0061 0.295 TRAF3, UFD1
GO:0001750 CC photoreceptor outer segment 23 2 0.123 0.0078 0.295 GNB1, IFT140
GO:0001655 BP urogenital system development 146 4 0.783 0.0084 0.295 ESR1, FOXF1, ID3, IFT140
GO:0000932 CC P-body 27 2 0.145 0.0104 0.295 AGO3, SYNE1
GO:0030522 BP intracellular receptor signaling pathway 85 3 0.456 0.0116 0.295 ESR1, ESRRG, UFD1

Meta Analysis With Previous Paper

The results in the Turretfield population from the paper Resistance to RHD virus in wild Australia rabbits: Comparison of susceptible and resistant individuals using a genomewide approach were then compared to these results. In this previous analysis, results were generated by using Fisher’s Exact Test to test genotypes for association with RHDV susceptibility or resistance.

tfGR <- url("https://raw.githubusercontent.com/hdetering/orycun/master/results/filteredSNPScores.csv") %>% 
  read_csv() %>%
  distinct(min_snp, .keep_all = TRUE) %>%
  dplyr::select(chromosome = chr_name, pos, snpID = min_snp, p, FDR) %>%
  mutate(chromosome = gsub("^chr", "", chromosome)) %>%
  makeGRangesFromDataFrame(
    keep.extra.columns = TRUE,
    ignore.strand = TRUE, 
    seqinfo = seqinfo(ensGenes),
    seqnames.field = "chromosome", 
    start.field = "pos", 
    end.field = "pos") %>%
  sort()

Of the 9,229 unique SNPs contained in the previous dataset, only 1,431 were also identified in the current dataset. A meta-analysis was performed on these SNPs using the results from both separate analyses which used Fisher’s Exact Test on genotypes. In this approach, Fisher’s method was used to combine p-values, and resultant p-values were drawn from a \(chi^2\) distribution with 4 degrees of freedom (i.e. \(2k\)).

combinedResults <- tfGR %>%
  find_overlaps(allSnpsGR) %>%
  mcols() %>%
  as.data.frame() %>%
  as_tibble() %>%
  dplyr::select(snpID = snpID.y, tf.p = p) %>%
  distinct(snpID, .keep_all = TRUE) %>%
  left_join(
    genotypeResults %>%
      dplyr::select(snpID, chromosome, position, geno.p = p)) %>% 
  mutate(
    X = -2*(log(tf.p) + log(geno.p)), 
    p = pchisq(X, 4, lower.tail = FALSE), 
    adjP = p.adjust(p, "bonferroni"),
    FDR = p.adjust(p, "fdr")) %>% 
  arrange(p)

Using Bonferroni’s adjustment to control the FWER at 0.05 gave 2 candidate SNPs, whilst using the more relaxed criterion of an FDR < 0.05 gave 8 candidate SNPs.

Significant SNPs from the meta-analysis using the previous published results from a separate population.
snpID Chromosome Position p adjP FDR
644062_1 3 56309981 1.668e-06 0.001615 0.001615
3378006_1 GL018704 2271135 2.355e-05 0.02279 0.0114
3941015_20 GL018847 216450 5.316e-05 0.05146 0.01715
213788_4 1 146274492 8.646e-05 0.0837 0.02092
3108287_38 20 10514947 0.0001406 0.1361 0.02722
3156991_0 21 6185579 0.0001792 0.1735 0.02892
3904662_77 GL018828 225354 0.0003003 0.2907 0.04014
3970912_82 GL018864 1306 0.0003317 0.3211 0.04014
Genes within 40kb of SNPs identified in the meta-analysis.
snpID Chromosome BP ID Name Description
213788_4 1 146274492 ENSOCUG00000000568 HBB2 hemoglobin, beta
213788_4 1 146274492 ENSOCUG00000021180 HBG2 hemoglobin, gamma G
213788_4 1 146274492 ENSOCUG00000027533 HBG1 hemoglobin, gamma A
213788_4 1 146274492 ENSOCUG00000026912 OR51I1 olfactory receptor family 51 subfamily I member 1
213788_4 1 146274492 ENSOCUG00000025404 OLFR64_1 MOR 5’beta3
644062_1 3 56309981 ENSOCUG00000004827 SFXN1 sideroflexin 1
3108287_38 20 10514947 ENSOCUG00000008859 KCNH5 potassium voltage-gated channel subfamily H member 5
3156991_0 21 6185579 ENSOCUG00000005961 KDM2B lysine demethylase 2B
3156991_0 21 6185579 ENSOCUG00000005957 RNF34 ring finger protein 34
3378006_1 GL018704 2271135 ENSOCUG00000011042 COL8A2 collagen type VIII alpha 2 chain
3378006_1 GL018704 2271135 ENSOCUG00000011038 ADPRHL2 ADP-ribosylhydrolase like 2
3378006_1 GL018704 2271135 ENSOCUG00000029722 TEKT2 tektin 2
3378006_1 GL018704 2271135 ENSOCUG00000006255 AGO3 argonaute RISC catalytic component 3
3904662_77 GL018828 225354 ENSOCUG00000027620 PAQR4 progestin and adipoQ receptor family member 4
3904662_77 GL018828 225354 ENSOCUG00000017181 PKMYT1 protein kinase, membrane associated tyrosine/threonine 1
3904662_77 GL018828 225354 ENSOCUG00000027985 TNFRSF12A TNF receptor superfamily member 12A
3904662_77 GL018828 225354 ENSOCUG00000017197 CLDN6 claudin 6
3941015_20 GL018847 216450 ENSOCUG00000012298 ODF2 outer dense fiber of sperm tails 2
3941015_20 GL018847 216450 ENSOCUG00000008855 CERCAM cerebral endothelial cell adhesion molecule
3970912_82 GL018864 1306 ENSOCUG00000021601 SCARB1 scavenger receptor class B member 1
Comparison of p-values for SNPs detected in both analyses, with thos considered as significant after meta-analysis shown in red.

Comparison of p-values for SNPs detected in both analyses, with thos considered as significant after meta-analysis shown in red.

Allele distributions in the 1996 and 2012 populations for SNPs identified under the meta-analysis.

Allele distributions in the 1996 and 2012 populations for SNPs identified under the meta-analysis.

GO Enrichment

This list was also used to test for GO enrichment.

sigSnpGR <- allSnpsGR %>%
  subset(snpID %in% filter(combinedResults, FDR < 0.05)$snpID)
nGenes <- allSnpsGR %>%
  subset(snpID %in% combinedResults$snpID) %>%
  unique() %>%
  length()
nSigGenes <- sigSnpGR$gene_id %>%
  unique() %>%
  length()
combinedGoRes <- bm %>%
  dplyr::filter(
    gene_id %in% subset(allSnpsGR, snpID %in% combinedResults$snpID)$gene_id) %>%
  mutate(inSig = gene_id %in% sigSnpGR$gene_id) %>%
  group_by(go_id) %>%
  summarise(N = dplyr::n(), 
            nSig = sum(inSig)) %>%
  left_join(goSummaries) %>%
  filter(nSig > 0, shortest_path > 2) %>%
  distinct(go_id, .keep_all = TRUE) %>%
  split(f = .$go_id) %>%
  mclapply(function(x){
    ft <- matrix(
      c(nSigGenes - x$nSig, x$nSig, nGenes - nSigGenes - x$N, x$N),
      nrow = 2, byrow = TRUE) %>%
      fisher.test()
    mutate(x, p = ft$p.value)
  }, mc.cores = mc) %>%
  bind_rows() %>%
  mutate(Expected = nSigGenes * N / nGenes,
         FDR = p.adjust(p, "fdr"),
         adjP = p.adjust(p, "bonferroni"),
         term = Term(go_id),
         ontology = Ontology(go_id)) %>%
  dplyr::select(go_id, term, ontology, shortest_path, N, nSig, Expected, p, adjP, FDR) %>%
  arrange(p)

Inspection of the results revealed an enrichment for terms related to haemoglobin, however, as these genes are co-located this was revealed to be due to a single SNP with multiple haemoglobin genes within 40kb and can be disregarded.

Most highly ranked GO terms for combined results. Where multiple terms resulted from the identical combination of genes, only the most highly ranked terms is shown. Overall, this corresponds to an FDR of 2.50%
ID Ontology Term Total Sig nSnps Expected p FDR Genes
GO:0015669 BP gas transport 4 3 1 0.07 2e-04 0.025 HBB2, HBG1, HBG2
GO:0019825 MF oxygen binding 4 3 1 0.07 2e-04 0.025 HBB2, HBG1, HBG2
GO:0015893 BP drug transport 14 4 2 0.245 2e-04 0.025 HBB2, HBG1, HBG2, SFXN1

Session Information

R version 3.6.1 (2019-07-05)

Platform: x86_64-pc-linux-gnu (64-bit)

locale: LC_CTYPE=en_AU.UTF-8, LC_NUMERIC=C, LC_TIME=en_AU.UTF-8, LC_COLLATE=en_AU.UTF-8, LC_MONETARY=en_AU.UTF-8, LC_MESSAGES=en_AU.UTF-8, LC_PAPER=en_AU.UTF-8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_AU.UTF-8 and LC_IDENTIFICATION=C

attached base packages: stats4, parallel, grid, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: GO.db(v.3.8.2), AnnotationDbi(v.1.46.1), Biobase(v.2.44.0), rtracklayer(v.1.44.4), plyranges(v.1.4.3), GenomicRanges(v.1.36.1), GenomeInfoDb(v.1.20.0), IRanges(v.2.18.2), S4Vectors(v.0.22.1), BiocGenerics(v.0.30.0), ggrepel(v.0.8.1), qqman(v.0.1.4), ggsn(v.0.5.0), rgdal(v.1.4-4), ggmap(v.3.0.0), sp(v.1.3-1), readxl(v.1.3.1), forcats(v.0.4.0), stringr(v.1.4.0), dplyr(v.0.8.3), purrr(v.0.3.2), readr(v.1.3.1), tidyr(v.0.8.3), tibble(v.2.1.3), ggplot2(v.3.2.1), tidyverse(v.1.2.1), magrittr(v.1.5), scales(v.1.0.0), pander(v.0.6.3), SNPRelate(v.1.18.1), SeqArray(v.1.24.2) and gdsfmt(v.1.20.0)

loaded via a namespace (and not attached): colorspace(v.1.4-1), rjson(v.0.2.20), class(v.7.3-15), futile.logger(v.1.4.3), XVector(v.0.24.0), rstudioapi(v.0.10), bit64(v.0.9-7), lubridate(v.1.7.4), xml2(v.1.2.2), knitr(v.1.24), zeallot(v.0.1.0), jsonlite(v.1.6), Rsamtools(v.2.0.0), broom(v.0.5.2), png(v.0.1-7), compiler(v.3.6.1), httr(v.1.4.1), backports(v.1.1.4), assertthat(v.0.2.1), Matrix(v.1.2-17), lazyeval(v.0.2.2), cli(v.1.1.0), formatR(v.1.7), htmltools(v.0.3.6), tools(v.3.6.1), gtable(v.0.3.0), glue(v.1.3.1), GenomeInfoDbData(v.1.2.1), Rcpp(v.1.0.2), cellranger(v.1.1.0), vctrs(v.0.2.0), Biostrings(v.2.52.0), nlme(v.3.1-141), xfun(v.0.9), rvest(v.0.3.4), XML(v.3.98-1.20), zlibbioc(v.1.30.0), hms(v.0.5.1), SummarizedExperiment(v.1.14.1), lambda.r(v.1.2.3), yaml(v.2.2.0), memoise(v.1.1.0), calibrate(v.1.7.2), stringi(v.1.4.3), RSQLite(v.2.1.2), highr(v.0.8), maptools(v.0.9-5), e1071(v.1.7-2), BiocParallel(v.1.18.1), RgoogleMaps(v.1.4.4), rlang(v.0.4.0), pkgconfig(v.2.0.2), bitops(v.1.0-6), matrixStats(v.0.55.0), evaluate(v.0.14), lattice(v.0.20-38), sf(v.0.7-7), labeling(v.0.3), GenomicAlignments(v.1.20.1), bit(v.1.1-14), tidyselect(v.0.2.5), plyr(v.1.8.4), R6(v.2.4.0), generics(v.0.0.2), DelayedArray(v.0.10.0), DBI(v.1.0.0), pillar(v.1.4.2), haven(v.2.1.1), foreign(v.0.8-72), withr(v.2.1.2), units(v.0.6-4), RCurl(v.1.95-4.12), modelr(v.0.1.5), crayon(v.1.3.4), futile.options(v.1.0.1), KernSmooth(v.2.23-15), rmarkdown(v.1.15), jpeg(v.0.1-8), blob(v.1.2.0), digest(v.0.6.20), classInt(v.0.4-1), VennDiagram(v.1.6.20) and munsell(v.0.5.0)

## [1] TRUE